% Calculate the the Error in G for the 1 sector model
%
% Jon Steinsson and Emi Nakamura, August 2006
%****************************************************

function GError = CalcGErrorCalvoPlus(G,F,F2,Q,a,rp,rad,alphaCalvo,...
    theta,Proba,ProbNgr)

agridsize = size(a,1);
rpgridsize = size(rp,1);
radgridsize = size(rad,1);
Ssize = agridsize*rpgridsize*radgridsize;

Q = reshape(Q,[Ssize 1]);
F = reshape(F,[Ssize 1]);
F2 = reshape(F2,[Ssize 1]);

% Get weights as a function of S_t/P_t-1
%********************************************************

QaInd = repmat((1:agridsize)',rpgridsize*radgridsize,1);

QradInd = repmat((1:radgridsize)',[1 agridsize rpgridsize]);
QradInd = permute(QradInd,[2 3 1]);
QradInd = reshape(QradInd,[Ssize 1]);

[dummy,QrpInd] = ismember(int32(F*10^6),int32(rp*10^6));
clear dummy
FInd = sub2ind([agridsize rpgridsize radgridsize],QaInd,QrpInd,QradInd);

[dummy,QrpInd] = ismember(int32(F2*10^6),int32(rp*10^6));
clear dummy
F2Ind = sub2ind([agridsize rpgridsize radgridsize],QaInd,QrpInd,QradInd);
clear QaInd QrpInd QradInd

OldInd = (1:Ssize)';

[FInd iFInd] = sort(FInd);
OldFInd = OldInd(iFInd);

[F2Ind iF2Ind] = sort(F2Ind);
OldF2Ind = OldInd(iF2Ind);

Q_old = Q;
% Apply F
Q = accumarray(FInd,Q_old(OldFInd));
Q = [Q; zeros(size(Q_old,1)-size(Q,1),1)];
Q2 = accumarray(F2Ind,Q_old(OldF2Ind));
Q2 = [Q2; zeros(size(Q_old,1)-size(Q2,1),1)];
Q = alphaCalvo*Q + (1-alphaCalvo)*Q2;

% Apply Proba
Q = ((reshape(Q,agridsize,rpgridsize*radgridsize)')*Proba)';

% Apply ProbNgr
Q = reshape(Q,[agridsize*rpgridsize radgridsize])*ProbNgr;
Q = reshape(Q,[Ssize 1]);

% Get the prices
%**********************************************************

GInd = CalculateGInd(a,rp,rad,G);

F = F(GInd);
F2 = F2(GInd);

% Create Price Index
%***************************************************************

F = exp(F).^(1-theta);
F2 = exp(F2).^(1-theta);

F = alphaCalvo*F + (1-alphaCalvo)*F2;

totWeight = sum(reshape(Q,[agridsize*rpgridsize radgridsize]))';
F = reshape(F.*Q,[agridsize*rpgridsize radgridsize]);

% Deal with zero mass
W_zero = (totWeight < 10^(-10));
totWeight = totWeight + W_zero;

sumF = sum(F)';
sumF_zero = (sumF < 10^(-10));
sumF = sumF + sumF_zero;

GError = log((sumF./totWeight).^(1/(1-theta)));
GError = (1-sumF_zero).*(1-W_zero).*GError ...
    + W_zero.*zeros(radgridsize,1) ...
    + (1-W_zero).*sumF_zero.*zeros(radgridsize,1);